/**
 * visTopoClustering.h
 * ubuntu8.04,g++4.3.2
 * 18/3/2009,VISUALIZATION_CGCAD_THSS_THU_BJ
 * Author : Xinlai,Lu
 *
 * Potential bugs:
 * The points on integration curves used to build the kd-tree may be not closed, ie. 
 * they may not reach the boundary and leave an opening, which may lead to that the
 * searching of sources in a domain access another different domain.
 *
 * Modification History:
 *
 */

#ifndef VIS_TOPO_CLUSTERING
#define VIS_TOPO_CLUSTERING

#include "ucdinputfile.h"
#include "visKDTree.h"
#include "visFlowTopology.h"
#include "visStreamlineClustering.h"
#include "visKMeans.h"
#include "visKMedoids.h"
#include "visHierarchicalClustering.h"

#include <ctime>
#include <algorithm>
#include <fstream>

#define SCANLINE_FILL

namespace visualization
{
//typedef double BOUNDARIES;
// One saddle is related to four domains(quadrants).
const size_t QUADRANTS = 4;
typedef vector<Point> SOURCES;
// Sources of four domains of a saddle.
typedef struct SourcesAroundSaddle
{
	SOURCES sources[QUADRANTS];
//	size_t centersNum[QUADRANTS];
}SourcesAroundSaddle;
// To save the clustering results.
typedef struct ClusterCenters
{
	vector<Point> *center;
	size_t cntNum;
}ClusterCenters;

class visTopoClustering : public visFlowTopology
{
public:
	visTopoClustering(Model *mdl) : visFlowTopology(mdl), m_ClusteringMethod(mdl),
																	m_SourcesInDomains(NULL), intgCurvesAndBoundaries(NULL),
																  m_ClusterCenters(NULL)
	{
		//m_FinalSourcesNum = 0;
		m_SelfIncrementCenters = 3;
		m_CrossIncrementCenters = 2;
#ifdef SCANLINE_FILL
		size_t dataPoint = mdl -> getPntNumberOfModel();
		mdl->getBound(xmin, xmax, ymin, ymax, zmin, zmax);
		dataPoint = size_t(sqrt((double)dataPoint));
		xres = yres = size_t((dataPoint + 1) * 2.5);
		m_xintv = (xmax - xmin)/xres;
		m_yintv = (ymax - ymin)/yres;
		m_pointNum = xres * yres;	
#endif
	}
	visTopoClustering() : visFlowTopology(), m_ClusteringMethod(),
												m_SourcesInDomains(NULL), intgCurvesAndBoundaries(NULL), 
												m_ClusterCenters(NULL)
	{
		//m_FinalSourcesNum = 0;
		m_SelfIncrementCenters = 3;
		m_CrossIncrementCenters = 2;
	}

	~visTopoClustering()
	{
		delete intgCurvesAndBoundaries;
		delete m_Mdl;
//		delete []m_ClusterCenters;
//		delete []m_SourcesInDomains;
		cout << endl << "REMEMBER TO DELETE THE MEMORY FOR : " << endl << 
						 				"-- m_SourcesInDomains" << endl <<
										"-- m_ClusterCenters" << endl << endl;
	}

	void setFileName(string file)
	{ m_FileName = file; }

	void parseParameter(string file);

	/**
	 *
	 */
	void preProcessSources();

	/**
	 * 
	 *
	 */
	void clusterDomain();

	/**
	 *
	 */
	void postProcess();
	void removeDuplicate();
	void makeArrow();

	/**
	 *
	 */
	void getTopoStructure(intgResults *&forward, intgResults *&backward,
												vector<saddlePoint> &saddles)
	{	
		forward = m_Forward;	
		backward = m_Backward; 
		saddles = getSaddles(m_Saddles);
	}

	void getBoundaries(double boundaries[6])
	{
		boundaries[0] = m_XMin;
		boundaries[1] = m_XMax;
		boundaries[2] = m_YMin;
		boundaries[3] = m_YMax;
		boundaries[4] = m_ZMin;
		boundaries[5] = m_ZMax;
	}

	SourcesAroundSaddle *getSourcesInDomain()
	{	return m_SourcesInDomains; }
	double m_XMin, m_XMax, m_YMin, m_YMax, m_ZMin, m_ZMax;

	ClusterCenters *getClusterCenters() {	return m_ClusterCenters; }
	size_t getCenterNumber(){ return m_CenterNum;}


	vector< vector<Point> > m_Streamlines;
	size_t m_CenterNum;
	size_t m_Iterations;
	size_t m_UseThreshold;
	size_t m_UseRemoveDuplicate;
	Point *m_Arrows;
private:
	/**
	 * The function determines the topology sturcture which will be used in function
	 * roamDomain() to build kdtree.
	 */
	void topoStructure();
	/**
	 * The function builds the kdtree using points on the boundary of the data field 
	 * and the integration curves.
	 */
	void buildKDTree();
	/**
	 * The function scans the saddles and it's related domains one by one, roams to 
	 * determine the sources in each domain.
	 */
	void roamDomain();
	SOURCES sourcesInDomain(Point &sourceSeed);
	void sourcesInDomain(SOURCES &srcs, Point sourceSeed, char direction);
//	void clearQueue()
//	{	for(size_t i = 0; i < m_PntQueue.size(); i++)	m_PntQueue.pop(); }

private:
	UCDInputFile *m_File;
	Model *m_Mdl;
	//visKMeans m_ClusteringMethod;
	//visKMedoids m_ClusteringMethod;
	visHierarchicalClustering m_ClusteringMethod;
	SourcesAroundSaddle *m_SourcesInDomains;
	size_t m_Saddles;
	//size_t m_FinalSourcesNum;
	kdTree *intgCurvesAndBoundaries;
	intgResults *m_Forward, *m_Backward;
	string m_FileName;
	ClusterCenters *m_ClusterCenters;
	size_t m_Segments;
	size_t m_MaxTime;
	double m_Step;
	double m_Scale;
	size_t m_SourcesNum, m_Interval;
	double m_Threshold;
	double m_DuplicateThshold;
	double m_DistanceThshold, m_Diavation;
	double m_GridLen;
	size_t m_GridInterval;
	
	double m_ClusteringThreshold;
	double m_preClusteringThreshold;
	/*Parameters related to clustering*/
	size_t m_SelfIncrementCenters;
	size_t m_CrossIncrementCenters;


#ifdef SCANLINE_FILL
	double xmin, xmax, ymin, ymax, zmin, zmax;
	size_t xres, yres, zres;
	double m_xintv, m_yintv, m_zintv;
	size_t m_pointNum;

	/// Test whether a point has been used as a seed to generate streamline.
	bool pointPassed(bool *matrix, size_t i , size_t j, size_t k)
	{
		if(2 == dim)
		{
		  if( i > xres || j > yres)
	 	  {
	 			cout << "Out of Boundary at line " << __LINE__ << " of file " << __FILE__ << endl;
	 			return true;
	 		}
	 		return matrix[j * xres + i];
		}
		else
		{
		  if( i > xres || j > yres || k > zres)
	 	  {
	 			cout << "Out of Boundary at line " << __LINE__ << " of file " << __FILE__ << endl;
	 			return true;
	 		}
	 		return matrix[k * (xres * yres) + j * xres + i];	
	  }
	
	}
	bool pointPassed(bool *matrix, Point &pnt)
	{
		size_t i, j, k;
		if(2 == dim)
		{
			i = size_t((pnt.getX() - xmin)/m_xintv);
			j = size_t((pnt.getY() - ymin)/m_yintv);
		}
		else
		{
			i = size_t((pnt.getX() - xmin)/m_xintv);
			j = size_t((pnt.getY() - ymin)/m_yintv);
			k = size_t((pnt.getZ() - zmin)/m_zintv);
		}
		return pointPassed(matrix, i, j, k);
	}

	void pointSet(bool *matrix, size_t i , size_t j, size_t k)
	{
		if(2 == dim)
		{
		  if( i > xres || j > yres)
	 	  {
	 			cout << "Out of Boundary at line " << __LINE__ << " of file " << __FILE__ << endl;
	 			exit(-1);
	 		}
	 		matrix[j * xres + i] = true;
		}
		else
		{
		  if( i > xres || j > yres || k > zres)
	 	  {
	 			cout << "Out of Boundary at line " << __LINE__ << " of file " << __FILE__ << endl;
	 			exit(-1);
	 		}
	 		matrix[k * (xres * yres) + j * xres + i] = true;	
	  }
	}

	void linePassingPoints(bool *matrix, Point &p1, Point &p2)
	{
		if(2 == dim)
		{
			size_t xstart = size_t((p1.getX() - xmin)/m_xintv);
			size_t ystart = size_t((p1.getY() - ymin)/m_yintv);
			size_t xend = size_t((p2.getX() - xmin)/m_xintv);
			size_t yend = size_t((p2.getY() - ymin)/m_yintv);
	
			size_t startX = xstart < xend ? xstart : xend;
			size_t endX = xstart > xend ? xstart : xend;
			size_t startY = ystart < yend ? ystart : yend;
			size_t endY = ystart > yend ? ystart : yend;
			for(size_t i = startX; i <= endX; i++)
			{
				for(size_t j = startY; j <= endY; j++)
				{
					pointSet(matrix, i,j,0);
				}
			}
	
		}
		else
		{
			size_t xstart = size_t((p1.getX() - xmin)/m_xintv);
			size_t ystart = size_t((p1.getY() - ymin)/m_yintv);
			size_t zstart = size_t((p1.getZ() - zmin)/m_zintv);
	
			size_t xend = size_t((p2.getX() - xmin)/m_xintv);
			size_t yend = size_t((p2.getY() - ymin)/m_yintv);
			size_t zend = size_t((p2.getZ() - zmin)/m_zintv);
	
			size_t startX = xstart < xend ? xstart : xend;
			size_t endX = xstart > xend ? xstart : xend;
			size_t startY = ystart < yend ? ystart : yend;
			size_t endY = ystart > yend ? ystart : yend;
			size_t startZ = zstart < zend ? zstart : zend;
			size_t endZ = zstart > zend ? zstart : zend;
			for(size_t i = startX; i <= endX; i++)
			{
				for(size_t j = startY; j <= endY; j++)
				{
					for(size_t k = startZ; k <= endZ; k++)
					{
						pointSet(matrix, i,j,k);
					}
				}
			}
	
	  }
	}
#endif // SCANLINE_FILL

};
}
#endif
